#### Measuring support for liberal democracy
#### Code for revised draft of paper
#### Israel sample
#### Version 5 - February 2024

library(tidyverse)
library(dplyr)
library(psych)
library(lavaan)
library(car)
library(semTable)
library(psychTools)
library(tth)
library(haven)
library(RColorBrewer)
library(matrixStats)
library(pewmethods)


# read stata file
sd_dat = read_dta("MeasuringDemocracy_IsraelData.dta")

# items
labelled::var_label(sd_dat)

# examine responses
names(sd_dat)
summary(sd_dat)
sapply(sd_dat[, 6:22], table, useNA="ifany") # check if DKs and NAs are included and how they are coded
# response set runs from 1=str disag to 5 str agr, 6 = DK
attr(sd_dat$FREXP1, "labels")

# drop DECELEC2 for this sample as the translated question 
# does not refer to an important *nonpolitical authority*
sd_dat = subset(sd_dat, select = -DECELEC2)


## Weights

# recode demogs
sapply(sd_dat[, c("age_group", "educ_level", "Gender", "Arab")], 
       function(x) table(x) / dim(sd_dat)[1] ) 
sapply(sd_dat[, c("age_group", "educ_level", "Gender", "Arab")], 
       function(x) attr(x, "labels") )
sd_dat = sd_dat %>% mutate(across(c(age_group, educ_level, Gender, Arab), as.factor))

# population data
pop_margins = list(
  tibble(age_group = as.factor(1:6), Freq = c(15.4, 20.1, 18.9, 15.5, 12.4, 17.8)),
  tibble(educ_level = as.factor(1:4), Freq = c(18.6, 41.5, 7.1, 33.0)),
  tibble(Gender = as.factor(c(0, 1)), Freq = c(48.8, 51.2)),
  tibble(Arab = as.factor(c(0, 1)), Freq = c(78.9, 21.1))
)

# raking
survey_weights = rake_survey(
  .data = sd_dat,
  pop_margins = pop_margins,
  scale_to_n = TRUE
)

# trim
sd_dat$wt_new = survey_weights
sd_dat$wt_new = trim_weights(survey_weights)
summary(sd_dat$wt_new)


## Plot item distributions

## Plot the distributions as the items were asked and answered, i.e., don't recode DKs or reflect the response set
## Include the response categories in this order: strongly agree, agree, neither, disagree, strongly disagree, DK

# colour palettes 
lik5 = c(brewer.pal(5, "BrBG"), "#b4b1b1")
lik5_labs = c("StrAg", "SomeAg", "Neither", "SomeDis",
              "StrDis", "DK")

# recode missing values as NAs; change "miss" to the missing values indicator; delete if already coded as NA
sd_plot = as.data.frame(sapply(sd_dat[, 6:21], 
                               function(x) car::recode(x, "1=5; 2=4; 3=3; 4=2; 5=1; 6=6; else=NA"), 
                               simplify = TRUE))

# plot
pdf("supdem_distr_israel.pdf", height=8, width=6)
par(mfrow=c(5, 4), mar=c(3, 2, 0.5, 0.5), tcl=-0.2, cex=0.9, las=2, mgp=c(1.8, 0.9, 0))
for(i in 1:length(sd_plot)) {
  barplot(height = table(sd_plot[, names(sd_plot)[i]]) / dim(sd_plot)[1] * 100, names.arg=lik5_labs, 
          axes=FALSE, col=lik5, cex.names=0.7, mgp=c(1, 0.2, 0), ylim=c(0, 80))
  axis(side=2, labels=TRUE, cex.axis=0.7, mgp=c(1, 0.4, 0))
  text(x=1.5, y=75, names(sd_plot)[i], cex=0.8, adj=0)
}
dev.off()


## Recode to orient all items in a pro-democratic direction and DKs as "neither"; 

sd_dat_r = sd_dat[, 6:21]
dem_vnam = c("FREXP1", "FRASSC2", "UNISUFF2", "FRELECT1", "JUDCNSTR1", "LEGCNSTR2", "EQLAW2")
aut_vnam = c("FREXP2", "FRASSC1", "FRASSC3", "UNISUFF1", "DECELEC1", "FRELECT2", "JUDCNSTR2", "LEGCNSTR1", "EQLAW1")
sd_dat_r[, aut_vnam] = lapply(sd_dat[, aut_vnam], car::Recode, "1=5; 2=4; 3=3; 4=2; 5=1; 6=3; else=NA")
sd_dat_r[, dem_vnam] = lapply(sd_dat[, dem_vnam], car::Recode, "1=1; 2=2; 3=3; 4=4; 5=5; 6=3; else=NA")
sd_dat[, aut_vnam] = sd_dat_r[, aut_vnam]
sd_dat[, dem_vnam] = sd_dat_r[, dem_vnam]


## Reliability and dimensionality

# alpha
sd_cor = cor(sd_dat_r, use="pair")
psych::alpha(sd_cor)
alph_out <- psych::alpha(sd_cor)
write.csv(alph_out[[1]], file="sd_alpha.csv")

# eigenvalues of principal components
eigen(sd_cor)$values
write.csv(eigen(sd_cor)$values, file="sd_eigen.csv")

# 1-factor EFA
psych::fa(sd_dat_r, nfactors=1)
efa1 = psych::fa(sd_dat_r, nfactors=1)
efa_html = tth(fa2latex(efa1))
writeLines(efa_html, "sd_efa1.html")

# ordinal
sd_pcor = polychoric(sd_dat_r)$rho
write.csv(sd_pcor, "sd_pcormat.csv", row.names=TRUE)
psych::fa(sd_pcor, nfactors=1, n.obs=dim(sd_dat)[1])
efa_ord = psych::fa(sd_pcor, nfactors=1, n.obs=dim(sd_dat)[1])
efa_ord_html = tth(fa2latex(efa_ord))
writeLines(efa_ord_html, "sd_efa1_ord.html")


## CFA models

# liberal democracy factor with orthogonal methods factor

cfa_mod_1 = '
SupLD =~ FREXP1 + FREXP2 + FRASSC1 + FRASSC2 + FRASSC3 + UNISUFF1 + UNISUFF2 
          + DECELEC1 + FRELECT1 + FRELECT2 + JUDCNSTR1 + JUDCNSTR2 
          + LEGCNSTR1 + LEGCNSTR2 + EQLAW1 + EQLAW2
PosVal =~ FREXP1 + FRASSC2 + UNISUFF2 + FRELECT1 + JUDCNSTR1 + LEGCNSTR2 + EQLAW2
'
sd_cfa_1_std = cfa(cfa_mod_1, data=sd_dat, estimator="MLR", orthogonal=TRUE, std.lv=TRUE, 
                   sampling.weights="wt_new")
semTable(sd_cfa_1_std, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="cfa_1fac_std_table_w.html")
sd_cfa_1_std_fit = fitMeasures(sd_cfa_1_std, output = "matrix",
                               fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_1_std_fit, file = "cfa_1fac_std_fit_w.csv", row.names = TRUE)

# liberal democracy factor with orthogonal methods factor, ordinal

sd_cfa_1_std = cfa(cfa_mod_1, data=sd_dat, estimator="WLSMV", orthogonal=TRUE, ordered=TRUE, std.lv=TRUE)
summary(sd_cfa_1_std, fit.measures=TRUE)
semTable(sd_cfa_1_std, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"), 
         type="html", file="cfa_1fac_ord_std_table.html")
sd_cfa_1_std_fit = fitMeasures(sd_cfa_1_std,  output = "matrix",
                               fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_1_std_fit, file = "cfa_1fac_std_ord_fit.csv", row.names = TRUE)

# electoral democracy and rule of laws factors with orthogonal methods factor

cfa_mod_2 = '
SupED =~ FREXP1 + FREXP2 + FRASSC1 + FRASSC2 + FRASSC3 + UNISUFF1 + UNISUFF2 
          + DECELEC1 + FRELECT1 + FRELECT2
SupRL =~  JUDCNSTR1 + JUDCNSTR2 + LEGCNSTR1 + LEGCNSTR2 + EQLAW1 + EQLAW2
PosVal =~ FREXP1 + FRASSC2 + UNISUFF2 + FRELECT1 + JUDCNSTR1 + LEGCNSTR2 + EQLAW2
PosVal ~~ 0*SupED
PosVal ~~ 0*SupRL
'
sd_cfa_2_std = cfa(cfa_mod_2, data=sd_dat, estimator="MLR", std.lv=TRUE, sampling.weights="wt_new")
semTable(sd_cfa_2_std, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="cfa_2fac_std_table_w.html")
sd_cfa_2_std_fit = fitMeasures(sd_cfa_2_std, output = "matrix",
                               fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_2_std_fit, file = "cfa_2fac_std_fit_w.csv", row.names = TRUE)

# electoral democracy and rule of laws factors with orthogonal methods factor

sd_cfa_2_std = cfa(cfa_mod_2, data=sd_dat, estimator="WLSMV", ordered=TRUE, 
                   std.lv=TRUE)
summary(sd_cfa_2_std, fit.measures=TRUE)
semTable(sd_cfa_2_std, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"), 
         type="html", file="cfa_2fac_ord_std_table.html")
sd_cfa_2_std_fit = fitMeasures(sd_cfa_2_std, output = "matrix",
                               fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_2_std_fit, file = "cfa_2fac_std_ord_fit.csv", row.names = TRUE)

# lr test - 1 v 2-fac models

lavTestLRT(sd_cfa_1_std, sd_cfa_2_std, method = "satorra.bentler.2010")
write.csv(lavTestLRT(sd_cfa_2_std, sd_cfa_1_std, method = "satorra.bentler.2010"), 
          file = "lrtest_cfas.csv", row.names = TRUE)


## Trimmed 7-item scale

sd_dat_7 = sd_dat_r[, c("FREXP2", "FRASSC1", "UNISUFF1", "FRELECT2", "JUDCNSTR2", "LEGCNSTR1", "EQLAW1")]

# 1-factor CFA

cfa_mod_7i = 'SupLD =~ FREXP2 + FRASSC1 + UNISUFF1 + FRELECT2 + JUDCNSTR2 + LEGCNSTR1 + EQLAW1'

sd_cfa_7i_std = cfa(cfa_mod_7i, data=sd_dat, estimator="MLR", std.lv=TRUE, 
                    sampling.weights="wt_new")
summary(sd_cfa_7i_std, fit.measures=TRUE)
semTable(sd_cfa_7i_std, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="cfa_7it_std_table_w.html")
sd_cfa_7i_fit = fitMeasures(sd_cfa_7i_std, output = "matrix",
                            fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_7i_fit, file = "cfa_7it_std_fit_w.csv", row.names = TRUE)

# 1-factor CFA, ordinal

sd_cfa_7i_ord = cfa(cfa_mod_7i, data=sd_dat, estimator="WLSMV", ordered=TRUE, std.lv=TRUE)
summary(sd_cfa_7i_ord, fit.measures=TRUE)
semTable(sd_cfa_7i_ord, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="cfa_7it_ord_std_table.html")
sd_cfa_7i_fit = fitMeasures(sd_cfa_7i_ord, output = "matrix",
                            fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(sd_cfa_7i_fit, file = "cfa_7it_ord_std_fit.csv", row.names = TRUE)

# eigenvalues of 7-item scale
sd7_cor = cor(sd_dat_7, use="pair")
eigen(sd7_cor)$values
write.csv(eigen(sd7_cor)$values, file="sd7_eigen.csv")

# reliability of 7-item scale
psych::alpha(sd7_cor)
alph7_out <- psych::alpha(sd7_cor)
write.csv(alph7_out[[1]], file="sd7_alpha.csv")

# ordinal reliability
sd7_pcor = polychoric(sd_dat_7)$rho
alph7_ord_out <- psych::alpha(sd7_pcor)
write.csv(alph7_ord_out[[1]], file="sd7_ord_alpha.csv")

# create additive 7-item scale
sd_dat$SUPDEM_7IT = rowMeans(sd_dat_7)


## Correlations with criterion variables

# create and save mixed correlation matrix

# linz question: recode such that support democ = 3, support auth = 1. 
table(sd_dat$Q5_Linzian, useNA="ifany")
attr(sd_dat$Q5_Linzian, "labels")
sd_dat$LINZ_SUPDEM = car::recode(sd_dat$Q5_Linzian, "1=3; 3=2; 2=1; else=NA")

# churchill question; recode such that support is high, oppose is low
table(sd_dat$Q6_Churchillian, useNA="ifany")
attr(sd_dat$Q6_Churchillian, "labels")
sd_dat$CHURCH_SUPDEM = car::recode(sd_dat$Q6_Churchillian, "8=NA")

# strong leaders question; recode such that support is high, oppose is low
table(sd_dat$Q7_WVS_1_Leader, useNA="ifany")
attr(sd_dat$Q7_WVS_1_Leader, "labels")
sd_dat$STRONG_LEAD = car::recode(sd_dat$Q7_WVS_1_Leader, "8=NA")

# importance of democracy question; recode such that support is high, oppose is low
table(sd_dat$Q8_ImportanceDemocracy, useNA="ifany")
attr(sd_dat$Q8_ImportanceDemocracy, "labels")
sd_dat$IMPORT_DEMOC = car::recode(sd_dat$Q8_ImportanceDemocracy, "11=NA")

# Satisfaction with democracy - SATIS_DEM: recode such that satisfaction is high
table(sd_dat$Q9_SatisfactionDemocracy, useNA="ifany")
attr(sd_dat$Q9_SatisfactionDemocracy, "labels")
sd_dat$SATIS_DEM = car::recode(sd_dat$Q9_SatisfactionDemocracy, "11=NA")

# Left-right / lib-cons ideology - LR_IDEOL or LC_IDEOL: recode such that right / conservative is high
table(sd_dat$Q16_IDEO, useNA="ifany")
attr(sd_dat$Q16_IDEO, "labels")
sd_dat$LR_IDEOL = car::recode(sd_dat$Q16_IDEO, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1; else=NA")

# Populist attitudes scale - POP_ATT
table(sd_dat$Q10_Populist_1_MK, useNA="ifany")
attr(sd_dat$Q10_Populist_1_MK, "labels")
sd_dat$pop1 = car::recode(sd_dat$Q10_Populist_1_MK, "6=3")
sd_dat$pop2 = car::recode(sd_dat$Q10_Populist_2_Public, "6=3")
sd_dat$pop3 = car::recode(sd_dat$Q10_Populist_3_Difference, "6=3")
sd_dat$pop4 = car::recode(sd_dat$Q10_Populist_4_Representation, "6=3")
sd_dat$pop5 = car::recode(sd_dat$Q10_Populist_5_Talking, "6=3")
sd_dat$pop6 = car::recode(sd_dat$Q10_Populist_6_Compromise, "6=3")
sd_dat$POP_ATT = rowMeans(sd_dat[, c("pop1", "pop2", "pop3", "pop4", "pop5")])

cor_items = c("SUPDEM_7IT", "LINZ_SUPDEM", "CHURCH_SUPDEM", "STRONG_LEAD", "IMPORT_DEMOC", 
              "SATIS_DEM", "LR_IDEOL", "POP_ATT")
crit_cor = mixedCor(sd_dat[, cor_items])$rho
write.csv(crit_cor, file="sd7_crit_cormat.csv")
